TP 2 Clasificacion

TP 2 Clasificacion. Tomas Palazzo y Axel Fridman

Cargamos las librerias

library("ggplot2")                  
library("GGally")

df = read.csv("lluviaAus.csv")

Limpieza y pre procesamiento de datos:

df$RainToday = as.factor(df$RainToday) # paso variables categoricas como factor
df$RainTomorrow = as.factor(df$RainTomorrow)

df$X <- NULL # borro columna X ya que sospechamos que no representa nada sino que es algun tipo de "id" que quedo grabado en el dataframe y no tiene influencia en la observacion. 

Chequeamos que cada columna sea del tipo correcto

str(df) 

'data.frame':   1000 obs. of  12 variables:
 $ MinTemp     : num  17.8 9 7.8 6.5 9 17.4 18.6 18.9 16.4 8.7 ...
 $ MaxTemp     : num  25.2 16 12.2 17.5 22.6 33.4 32.6 35.5 22.9 22.1 ...
 $ Rainfall    : num  0 0.8 1.8 0 0 0 0 0 0 0 ...
 $ Evaporation : num  4 1.6 0.6 2 2.8 6.8 9 15.2 3.6 3 ...
 $ Sunshine    : num  6.4 7.4 5.4 9.7 9.5 10.5 11.8 12.2 10.4 8 ...
 $ WindSpeed3pm: int  13 26 22 7 11 20 15 31 13 15 ...
 $ Humidity3pm : int  66 53 84 47 42 18 38 13 57 39 ...
 $ Pressure3pm : num  1013 1013 1026 1025 1016 ...
 $ Cloud3pm    : int  7 7 6 2 5 1 1 1 2 5 ...
 $ Temp3pm     : num  24.4 14.8 10.5 17.2 21.6 32.1 29.8 34.9 21.8 21.5 ...
 $ RainToday   : Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
 $ RainTomorrow: Factor w/ 2 levels "No","Yes": 1 1 2 1 2 1 1 1 1 1 ...

Analisis exploratorio de datos y visualizaciones:

g = ggpairs(df, progress = FALSE, bins=10)+theme_bw()
g

Lo que observamos es que hay variables que estan muy correlacionadas linealmente, como la de temperatura maxima y temperatura a las 3pm. Y otras que parecen ser muy poco relacionadas como la humedad y la presion.

print(paste("Correlacion humedad y presion (baja) " , cor(df$Humidity3pm, df$Pressure3pm) ))

[1] "Correlacion humedad y presion (baja)  0.0456887832592021"

print(paste("Correlacion maxTemp y temp a 3pm (muy alta)" , cor(df$Temp3pm, df$MaxTemp) ))

[1] "Correlacion maxTemp y temp a 3pm (muy alta) 0.980580835987693"

Tambien notamos que aproximadamente 4/5 de las observaciones no llueve, tanto ese mismo dia como el siguiente.

table(df$RainToday)

 No Yes 
802 198 

table(df$RainTomorrow)

 No Yes 
794 206 

Pero no son para nada independientes si llueve hoy y si llueve mañana. Ya que, si solo tuviera de informacion estas 2 columnas, dado que llovio hoy la nueva probabilidad de que llueva mañana es aproximadamente 78/(120+78) =aprox 40% mucho mas que el 20% naive.

dfLluvias = df[c("RainToday", "RainTomorrow")]
table(dfLluvias)

         RainTomorrow
RainToday  No Yes
      No  674 128
      Yes 120  78

Ejercicio 2

ggplot(df, aes(x=Sunshine, y=Humidity3pm, color=RainTomorrow)) +
  geom_point() 

Observamos cosas muy relevantes, los dias que va a llover mañana, tienen mayor humedad y menos sol, y los dias que no llovera mañana tienen todos mucho mas sol y la humedad en promedio es mas baja. A su vez hay varios dias, en su mayoria dias que llovera mañana, cuyo nivel de sol es 0, lo cual genera esa columna en el lado izquierdo. Tambien vemos que si bien hay muchos dias que tienen mucho sol y poca humedad (los dias que no llovera mañana), no vemos casi ninguna observacion con poco sol y poca humedad, lo cual nos podria hablar de cierta relacion humedad - sol. Nosotros pensamos que la humedad tiene mayor capacidad predictiva (si se tomara una sola y no en conjunto ambas), ya que la variable sunshine esta mucho mas dispersa para los dias que No llovio al dia siguiente, lo analizaremos en dos graficos de densidad.

ps<-ggplot(df, aes(x=Sunshine, fill=RainTomorrow)) +
  geom_density(alpha=0.4) + labs(x= "Nivel de radiacion solar (Sunshine)",
       subtitle="Grafico densidad radiacion solar") + geom_vline(xintercept=7.8, size=0.5, color="red")
ps

ph<-ggplot(df, aes(x=Humidity3pm, fill=RainTomorrow)) +
  geom_density(alpha=0.4) + labs(x= "Nivel de humedad (Humidity3pm)",
       subtitle="Grafico densidad humedad") + geom_vline(xintercept=62.5, size=0.5, color="red")

ph

Despues de ver estos 2 graficos notamos que no es facil dar un punto de corte para diferenciar las 2 clases solamente tomando una variable. Ya que si tomasemos aproximadamente 7.8 de punto de corte para la radiacion solar o 62 para la humedad como punto de corte, de todas formas tendrias bastante error ya que las clases se solapan mucho mirandolas unidimensionalmente. Tomamos estos valores de referencia como para decir que ningun corte es bueno, estos ni siquiera tienen en cuenta la diferencia de proporcion de clases. En definitiva no nos casamos con ninguna variable.

Ejercicio 3. Los boxplot no son buenos graficos. Pueden ocultar demasiada informacion cuando hoy en dia tenemos la capacidad de procesarla.

p2<-ggplot(df, aes(y=Humidity3pm, x=RainTomorrow, fill=RainTomorrow)) +
  geom_boxplot() + labs(x= "Nivel de humedad (Humidity3pm)",
       subtitle="Boxplots humedad segun lluvia mañana")
p2

p3<-ggplot(df, aes(y=Sunshine, x=RainTomorrow, fill=RainTomorrow)) +
  geom_boxplot() + labs(x= "Nivel de radiacion solar (Sunshine)",
       subtitle="Boxplots radiacion solar segun lluvia mañana")
p3

Como ven las medias difieren para ambos casos, pero esa informacion ya la habiamos visto (ademas de muchas otras cosas que aca no) en los density plots. No hay outliers / valores atipicos.

Ejercicio 4 Para hacer las ventanas moviles voy a primero transformar el dataset en 1 y 0 a las categoricas, para poder luego tomar promedios.

df$RainToday = ifelse(df$RainToday=="Yes",1,0)
df$RainTomorrow = ifelse(df$RainTomorrow=="Yes",1,0)

promediosMoviles = function(datosX, datosY, alturaAestimar, h){
  df = data.frame(datosX,datosY)
  df2 = df[df$datosX > alturaAestimar-h & df$datosX < alturaAestimar+h ,]
  if(mean(df2$datosY)>1/2 ){
    return(1)
  }
  return(0)
}

promediosMoviles(df$Sunshine, df$RainTomorrow, 8, 0.1)

[1] 0

promediosMoviles(df$Humidity3pm, df$RainTomorrow, 85, 2)

[1] 1

Ejercicio 5 y 6. Nos creamos una funcion que nos genere todo el vector de predicciones Yhat.

generarColumnaPrediccionesPromediosMoviles = function(datosX, datosY , h){
  sexoPredicho = c()
  for (i in (1: length(datosX))){
    sexoPredicho[i] = promediosMoviles(datosX, datosY, (datosX[i]), h)
  }
  return(sexoPredicho)
}

(generarColumnaPrediccionesPromediosMoviles(df$Sunshine, df$RainTomorrow, 1))

   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0
  [45] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  [89] 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0
 [133] 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0
 [177] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
 [221] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [265] 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [309] 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
 [353] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [397] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
 [441] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 [485] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0
 [529] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0
 [573] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [617] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [661] 1 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
 [705] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0
 [749] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1
 [793] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [837] 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0
 [881] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
 [925] 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
 [969] 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Vamos a evaluarlo con el metodo de validacion cruzada de LOO (dejar uno afuera para entrenar y evaluarlo con ese).

leaveOneOut = function(datosX, datosY, h){
  error = 0
  for (i in (1: length(datosX))){
    predichoI = promediosMoviles(df$Sunshine[-i], df$RainTomorrow[-i], df$Sunshine[i], h)
    error = error + abs(predichoI - df$RainTomorrow[i])
  }
  return(error)
}

hPosibleHumedad = seq(0.2, 10, 0.8 )
hPosibleSunshine = seq(0.5, 4, 0.1 )

erroreshHumedad = c()
for (i in (1: length(hPosibleHumedad))){
    erroreshHumedad[i] = leaveOneOut(df$Humidity3pm, df$RainTomorrow, hPosibleHumedad[i])
}

erroreshSunshine = c()
for (i in (1: length(hPosibleSunshine))){
    erroreshSunshine[i] = leaveOneOut(df$Sunshine, df$RainTomorrow, hPosibleSunshine[i])
}

plot(hPosibleHumedad , erroreshHumedad, type = "l")

De aca vemos que la ventana optima para humedad es de exactamente 1.

plot(hPosibleSunshine , erroreshSunshine, type = "l")

Mientras que la ventana optima para la radiacion solar es de 0.8

Asi el error de clasificacion empirico despues de encontrar el h optimo para sunshine es de 179. Es decir le erra el 17.9% de las veces